function bdy3

%  solve bvp

clear *
clf

nx=200;
% get(gcf)
set(gcf,'Position', [1281 735 577 197]);

hold on

ep=0.1;
dy=1+0.9*1;
set(gca,'ytick',[1  2]);

%ep=0.01;
%dy=1+0.88*2;
%set(gca,'ytick',[1 2 3]);

% find composite solution
xx=linspace(0,1,nx+2);
r=(-1+sqrt(5))/2;
for ix=1:nx+2
	composite(ix)=exp(xx(ix))+exp(-xx(ix)/ep)+(1-exp(1))*exp(r*(xx(ix)-1)/ep);
end;

%  find numerical solution
x0 = 0.0;	y0 = 2.0;
x1 = 1.0;	y1 = 1.0;

% parameters for calculation
nx = 400;
error = 0.000001;
x=linspace(0,1,nx+2);
for ix=1:nx+2
	y(ix)=y0+(y1-y0)*(x(ix)-x0);
end;
dx=x(2)-x(1);
dxx = dx*dx;
err=1;
counter=0;
while err > error
	% calculate the coefficients of finite difference equation
	a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
	for j = 2:nx+1
		jj=j-1;
		z = (y(j+1) - y(j-1))/(2*dx);
		a(jj) = 2 + dxx*fy(x(j), y(j), z,ep);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z,ep);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z,ep);
	end;
	% Newton iteration
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:nx
		xl = a(j) - c(j)*u(j-1);
		v(j) = (v(j) - c(j)*v(j-1))/xl;
		u(j) = - (2 + c(j))/xl;
	end;
	vv = v(nx);
	y(nx+1) = y(nx+1) + vv;
	err = abs(vv);
	for jj = nx:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
	end;
	counter=counter+1;
end;

%  plot results
plot(x,y,'k','LineWidth',1)
hold on

plot(xx,composite,'--r','LineWidth',1)
box on
%axis([0.01 1 0.5 1.1])
%legend(' Numerical',' Composite',4);

legend1 = legend(' Numerical',' Composite');
set(legend1,'Position',[0.6783 0.3133 0.118 0.09836]);

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
text(0.1,dy,'\epsilon = 0.1','FontSize',16,'FontWeight','bold');
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off


function g=f(x,y,z,ep)
g=(-exp(x)-ep*x*z+y)/ep^2;

function g=fy(x,y,z,ep)
g=1/ep^2;

function g=fz(x,y,z,ep)
g=-x/ep;


